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Abstract 

Vortices in the simplest superconducting state of graphene contain very low energy excita- 
tions, whose existence is connected to an index theorem that applies strictly to an approximate 
form of the relevant Bogoliubov-deGennes equations. When Zeeman interactions are taken 
into account, the zero modes required by the index theorem are (slightly) displaced. Thus the 
vortices acquire internal structure; the resulting "modicules" obey nonabelian quantum statis- 
tics. 
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In this paper we will draw together several lines of thought. Electronic properties of the 
two-dimensional material graphene have attracted theoretical interest for many years JU, and 
of course recently [2]. Many of graphene's unusual properties derive from the fact that its con- 
duction and valence bands touch at two points, forming conical energy surfaces near those 
points. At neutral filling, the Fermi energy coincides with the apex of the cone. As we shall 
discuss momentarily, one can induce superconductivity in graphene, probably in several forms 
||3l . Even the simplest such superconducting state has been found to have unusual properties 
with respect to Andreev reflection J4[. Here we point out that vortices, or more generally mul- 
tivortices, in this state acquire interesting internal structure. This occurs because each vortex 
supports a low-energy mode of the equation for electronic excitations, i.e. the Bogoliubov- 
deGennes (BdG) equations [5 ]. Indeed, an approximate form of the BdG equation maps, after 
appropriate identifications, to an equation of relativistic field theory that has been investigated 
by Jackiw and Rossi who discovered zero-energy solutions. E. Weinberg [7 ] subsequently 
demonstrated that the existence of these zero-energy solutions is connected to an index the- 
orem. The structured vortices resemble in some respects vortices in p + ip superconductors 
|[8l or in the Pfaffian quantum Hall state [9]. In particular, they support nonabelian statistics 
iflOl . a phenomenon of fundamental interest that has also played a central role in proposals for 
topological quantum computing ifTTTl . 

Induced Superconducivity in Graphene 

At neutral filling the density of states at the Fermi surface, which degenerates to two points, 
is very small, and in two dimensions fluctuations are important lfT2l . so the prospect for in- 
trinsic superconductivity in undoped graphene is problematic. If, however a graphene layer 
is put in contact with a superconducting substrate, then electron-electron interactions can in- 
duce anomalous (non-conserving) terms in the effective Hamiltonian, according to the general 
scheme 

n = -g^&W + h c _, _ -g^yi)^ + h.c. = -gA^-ip + h.c. (1) 

where \& is the total electron field, Ao is the bulk condensate, and tp is the electron field in 
the graphene layer; here only the anomalous terms have been retained. Upon diagonalizing the 
quadratic Hamiltonian for tp, a graphene condensate A = (iptp) is induced. 

That broad-brush sketch took no notice of the intricate internal structure of the graphene 
field tp and (possibly) of Ao and consequently A. The relevant, low-energy modes of tp are 
labeled by a 2-momentum k and three binary indices. The first binary index is the valley index, 
which specifies whether the mode arises from expansion around total momentum p = ±K + 
k; here ±K are the two momenta where the bands touch. The second binary index labels 
a pseudospin, arising from the non-trivial residual symmetry of the unit cell, which roughly 
speaking specifies whether the electron is on the A or B sublattice of the bipartite honeycomb 



lattice. This pseudospin is very significant dynamically, as it appears in the effective Dirac 
equation for these modes (see below). Finally, the third binary index labels ordinary spin. 

As in helium-3 lfT3l . or for that matter QCD [ 14], internal structure for the fermion field 
opens up many possibilities for the form of condensation. For definiteness, let us assume that 
the bulk condensate pairs electrons of opposite total 2-momentum. (This excludes LOFF-type 
bulk superconductivity |fT51 .) Then only intervalley pairing of the form 

A + ** ; _**(A;) = (V>+** (&)-!/>-**(-&)) (2) 

is induced; intravalley pairing requires Aq(±2K) ^ 0. (The possibility of intrinsic intravalley 
pairing has been discussed [16].) Here the three indices are the respective binary indices men- 
tioned previously and * is a wildcard. A general restriction arises from Fermi statistics, but it is 
very weak: components of the condensate which are overall antisymmetric under interchange 
of spatial, pseudospin, and spin must be symmetric in the valley index, and vice versa. 

With different bulk superconductors, one can imagine many exotic possibilities being re- 
alized, e.g. gapless or gapped p-wave or d-wave pairing or LOFF states lfl31 . The most con- 
ventional superconductors, however, are T-invariant, s-wave, and spin singlet. Assuming that 
the primary interaction in Eqn. (Q]) conserves spin (as is appropriate at least for light elements), 
then the induced superconductivity will likewise be s-wave spin singlet, and therefore sym- 
metric in space and antisymmetric in spin. 

The underlying Cq v symmetry of the honeycomb lattice on a homogeneous substrate im- 
poses a further restriction. With the valley index fixed, we have a residual C$ v . Consulting the 
character tables of the double group of C% v ifTTl . we find that the singlet piece of the tensor 
product of any non-trivial two-dimensional representation with itself lies in the symmetric part; 
and of course it combines pseudospin up and down. Overall antisymmetry then also requires 
symmetry in the valley indices; similarly, we already argued they should be different. 

The preceding discussion pointed to many byways worthy of further investigation. For our 
present purpose, however, the central conclusion is that conventional bulk superconductors will 
induce a very specific form of condensate, describable by a single complex-number field. We 
shall adopt that choice, implicit in 0, in what follows. 

Multivortices and Near-zero Modes 

The BdG equations from (4J, extended to include an electromagnetic vector potential, take 
the form 
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where we absorb the coupling constant into A, put h = Vf = 1, measure energies relative to 
the cone apices, and write H± = H± + H±, defining H± = —i(a x d x ± a y d y ) and H± = 
—q(a x A x ± a y Ay). Here the subscripts refer to the valley index and the internal indices, the 
internal indices are for pseudospin, u, v refer to particle and hole modes, respectively. In this 
approximation ordinary spin is taken to be dynamically inert, and does not appear. Note that 
the condensate mixes spin up particles with spin down holes, and vice versa. 

These equations decouple into two independent sets, one involving u_) and H + , the 
other those objects with the complementary indices. The equations for the second set can be 
related to those for the first by reflection about the x axis. Restricting to the first set, dropping 
the indices, and focusing on E = 0, we find the equations 

[-i(a x d x + aydy) - q(a x A x + a y A y )) u + A(r) v = (4) 
A*(r) u + [i(cr x d x + a y d y ) - q(a x A x + <r y A y )) v = (5) 

Now putting v = a y u* we find that the second equation reduces to the complex conjugate of 
the first, which reads 

[-i(a x d x + a y d y ) - q{a x A x + a y A y )] u + A(r) a y u* = (6) 

Equation © is precisely the equation for zero modes of the relativistic two-component Dirac 
equation studied in Refs. @|7). 

For later use, and to make our discussion self-contained, let us briefly review the solutions 
of Eqn. © for multivortices. With A = — egA(r), and using polar coordinates, we find for the 
upper component of u - call it a - the equation 

e ^l" + -4)° " + Aa* = (7) 

or r oO 

and for the lower component, b, the equation 

" e " l6 (|- " ~jL) b + V Ae ~ idb - A6* = (8) 
or r ot> 

Let us focus on the former. First, we can remove the vector potential term by rescaling a = 
exp(q Jq A(s)ds)a. For a vortex A vanishes at the origin; it also vanishes as r _1 as r — > do, 
so this redefinition will not affect the normalizability of the zero-modes below, which die 
exponentially at infinity. 

We suppose that A(r,9) = A(„)(r)e m6 ', with A( n )(r) — > r n l as r — > and A( n \(r) — > 
const, as r — > oo, as is appropriate to an ra-fold multivortex. (In this situation A(r) — > 
so the vector potential asymptotically "pulls in" a by a power of half the vorticity.) The def- 
inite partial wave solutions of Eqn. (O are of two types. One type involves a single angular 
dependence, a = f(r)e iW . The consistency condition for I is easily found to be 21 = n — 1. 
The second type involves two angular dependencies, a = f (r)e tW + g(r)e tm6 . The consistency 
condition for I and m is I + m = n — 1. 
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For the first type, we derive the radial equation 

f'- l -f + ^n)f* = (9) 

Assuming for simplicity that Ar n -\ is real, and positive at infinity, / can be taken real. At large 
r the middle term can be dropped, and the solution dies exponentially. At small r the last term 
can be dropped, and we see that solution is normalizable for I > — |, or, since I is integral, 
I > 0. For the second type, we derive the radial equations 

/' - l ~f + \n)9* =0 (10) 

TTl 

g'--g + A {n) f* = (11) 

Now there are both growing and dying modes at infinity. In order to assure that the dying mode 
matches onto a normalizable solution at the origin, we must require both l,m > 0. With the 
same assumptions on A( n ), there is a solution /o, #o with / an d g both real and in phase at 
infinity, and another i/o, —igo with / and g pure imaginary and out of phase. On the other 
hand, from the definition interchange of I and m is a trivial operation. Thus for positive n 
we find an n-real dimensional manifold of zero modes. When n is odd there is one solution 
of the first type, corresponding to I = and n — 1 of the second type, corresponding to 
n — 1 > I > 1. When n is even all solutions are of the second type, with n — 1 > I > 0. 

A very similar analysis applies to the b equation, Eqn. ([U). In that case, there are \n\ nor- 
malizable modes for n < 0, and none for n > 0. In all cases, the number of solutions of 
the a equations minus the number of solutions of the b equations equals n, the vorticity. This 
suggests the existence of an underlying index theorem, since index theorems generally express 
the difference between the number of zero modes of a differential operator V and the number 
of solutions of its adjoint - the index of V - in terms of topological data in the structure of 
V iTTSl . Index theorems can be very valuable in physics, because they reveal the presence of 
low-energy modes whose existence might otherwise be hard to anticipate, and insure that the 
existence of such low-energy modes is robust against many kinds of perturbations (or approxi- 
mations). Following E. Weinberg, if we write for the upper component a = h + ik, with h and 
k real, and similarly A = (pi + ifo, and combine h, k into a two-component spinor, then the 
Dirac equation Eqn. © takes the form 




while for the lower component, similarly decomposed, we find the adjoint equation. There is 
indeed an index theorem for V, essentially equating the index to the vorticity f7l . 

The equations involving (u-,v + ) and H- contribute another \n\ zero modes, after a par- 
allel analysis. Finally there is another overall doubling, when we restore the physical spin 
variable. 



So far our discussion has been based on the approximate BdG equation (0>, which does not 
include the Zeeman coupling of spin (as opposed to pseudospin) to the magnetic field. This is 
a small effect quantitatively, but it has significant qualitative and conceptual implications. The 
Zeeman coupling makes an additional diagonal contribution ±nBl to the matrix in ©, where 
B is the magnetic field, with the sign depending on spin. (Intuitively: since the background 
Cooper pairs are spin singlets, mixing with them does not affect the Zeeman energy.) As a 
result the former zero-energy states are shifted. In first-order perturbation theory in n, we have 
the shifts 

e± = ±k J0 ";T,r V ' ' (13) 



/ °° dr2nrB(r)\a(r) |2 
Jq 00 dr2Trr\a(r)\ 2 



Thus it is roughly proportional to the average magnetic field over the mode. For a crude es- 
timate of the "minimal" splitting, take the product of the bare ^-factor of an electron and the 
quantized minimal quantum fluxoid spread over a square micron, then e ~ 10~ 7 eV, i.e. 1 mK. 
Larger splittings might be obtained using in-plane magnetic fields, which couple to spin but do 
not frustrate the superconducting order parameter. 

It might seem paradoxical that a small perturbation can shift — and thus remove, as strict 
zero modes - zero modes whose existence was tied to topology. In the present context, how- 
ever, inclusion of the Zeeman term blocks passage from the amended Eqns. (|4j [5]) to Eqn. 
©, for which the index theorem applies. Alternatively, we could formally include a Zeeman 
coupling-like term directly in Eqn. ((6]). We would then still find zero modes to satisfy the in- 
dex theorem for the slightly perturbed V, but we could not use them to get zero modes of the 
amended Eqns. (01 [5). 

Comments 



1. It is convenient to have a single word to convey the concept "soliton that acquires internal 
structure due to existence of localized low-energy modes of quantum fields it interacts 
with". We propose the term modicule, pronounced mode-icule, in view of the resem- 
blance of such entities to emergent molecules. 

2. Experimental probes: General techniques for probing the internal structure of vortices 
and multivortices, notably including tunneling microscopy, were outlined in the pioneer- 
ing work of Virtanen and Salomaa fT9ll . As they emphasized, multivortices can be en- 
couraged to form at defects. The characteristic Zeeman splittings, predicted above, might 
be probed in absorption. (Note that the spatial wave functions for opposite spins are ac- 
curately matched.) 

3. Ragged multivortices and edges: Because the existence and number of zero modes is 
governed by an index theorem, and the relevant topological information is insensitive 
to the existence even of large holes where A vanishes, we will have near-zero modes 



associated with the edge of an annular region of graphene superconductor threaded by 
magnetic flux. Nor is symmetry required. 

4. Nonabelian statistics: For odd n, and in particular for the unit vortex n = ±1, we have 
unpaired, essentially real ("Majorana") modes. There are four of them, due to the spin 
and valley degeneracies. To a first approximation they are dynamically independent, 
because the spin couples feebly and intervalley scattering requires large momentum 
exchange ±2K. If we quantize them separately we find, following Ivanov ifTUl . non- 
abelian statistics: essentially, four copies of the Clifford algebra discovered by Nayak 
and Wilczek [20]. Whether and how this structure is modified in higher approximations 
is a problem worthy of much further investigation. Superconducting vortices can be ma- 
nipulated electrically, so there are good prospects for controlled interchanges of the mod- 
icules, enabling experiments to test whether such interchanges affect their structure in the 
predicted fashion. We will report some specific ideas along this line elsewhere CTTl . Note 
that the spin-dependent interactions allow one, in principle, to perform additional opera- 
tions on the states, such as (specifically) shifting relative phases in a manner controlled 
by magnetic fields. This additional control could be helpful in quantum information pro- 
cessing (22l. 

5. Comparing p + ip: An effective Dirac equation of similar structure appears in the the- 
ory of superconductors with bulk p x + ip y pairing [8]. The approximation underlying 
it, however, is quite different. Specifically, the momentum dependence is assumed to 
arise from the gap parameter, by expanding locally A(r, p) — > \A(r)\e l9 ( r \p x + ip y ) in 
the Bogoliubov-deGennes equation, where 6 is the phase of the order parameter. There 
are two potential problems with this approximation: it applies only if the ordinary mo- 
mentum dependence can be neglected, and it breaks down when [ A(r) | vanishes. The 
first assumption is dubious, when applied to the low-energy modes near a generic Fermi 
surface; and the second precludes clean application to vortices or edges. 
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